{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kinetic Mechanism Pathway Analyzer\n",
    "\n",
    "This analyzer is designed to provide a better user experience and improved analyzing capacity (e.g., large mechanisms) compared with Chemkin mechanism analyzer. This analyzer takes reaction ROP file (e.g. Chemkin ckcsv file), builds a reaction network, helps answer\n",
    "\n",
    "- if the mechanism has certain pathway (from start species to end species)\n",
    "\n",
    "- what is the biggest flux from/towards a certain species\n",
    "\n",
    "When using this script, please pay attention to the sections marked with **[user input]** which require user to specify certain input parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rmgpy.rmg.model import CoreEdgeReactionModel\n",
    "from rmgpy.chemkin import load_chemkin_file\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Prepare files\n",
    "\n",
    "Please add \n",
    "\n",
    "- annotated chemkin files (mechanism and species dictionary) and \n",
    "\n",
    "- chemkin simulation ckcsv file \n",
    "\n",
    "to a folder with name stored in `mech` variable in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mech = 'minimal'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mech_path = os.path.join('./', 'data', 'pathway_analysis', mech)\n",
    "chemkin_path= os.path.join(mech_path, 'chem.inp')\n",
    "dictionary_path = os.path.join(mech_path, 'species_dictionary.txt')\n",
    "ckcsv_path= os.path.join(mech_path, 'CKSoln.ckcsv')\n",
    "\n",
    "model = CoreEdgeReactionModel()\n",
    "model.core.species, model.core.reactions = load_chemkin_file(chemkin_path, dictionary_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate pairs for reactions that don't have flux pairs\n",
    "for rxn in model.core.reactions:\n",
    "    if not rxn.pairs:\n",
    "        rxn.generate_pairs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create reaction network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from rmgpy.tools.ckcsvparser import get_concentration_dict_from_ckcsv, get_rop_from_ckcsv, get_flux_graph_edges_dict\n",
    "from rmgpy.chemkin import get_species_identifier\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "from rmgpy.rmg.pdep import PDepReaction\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### step1: prepare a dict for graph edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_col_dict, spc_mf_dict = get_concentration_dict_from_ckcsv(ckcsv_path)\n",
    "first_col_dict_rop, spc_total_rop_dict, spc_rop_dict = get_rop_from_ckcsv(ckcsv_path)\n",
    "\n",
    "graph_edges_dict = get_flux_graph_edges_dict(spc_rop_dict, model.core.reactions)\n",
    "\n",
    "graph_edges_dict_simple = {}\n",
    "for pair in graph_edges_dict:\n",
    "    node1 = get_species_identifier(pair[0])\n",
    "    node2 = get_species_identifier(pair[1])\n",
    "    graph_edges_dict_simple[(node1, node2)] = graph_edges_dict[pair]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [user input] step2: pick a time to analyze pathways\n",
    "\n",
    "Please specify \n",
    "\n",
    "- the time point you want to investigate\n",
    "\n",
    "- the chemkin name of your initial species/reactant\n",
    "\n",
    "- the reaction temperature in K\n",
    "\n",
    "- the reaction pressure in Pa if the mechanism includes pressure-dependent reactions\n",
    "\n",
    "- show_integrated_flux: True if users want to use time-integrated fluxes (integrated up to the investigated time point) in mole for analysis. Flase if users want to use instantaneous fluxes in mole/cm3/s at the investigated time point for analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "time_investigated = 0.000001 # hour\n",
    "reactant_name = 'ethane(1)'\n",
    "T = 1350 # K\n",
    "P = 100000 # Pa\n",
    "show_integrated_flux = True \n",
    "##########################################\n",
    "\n",
    "\n",
    "timepoint_index = (np.abs(first_col_dict['Time_(sec)']-time_investigated*3600)).argmin()\n",
    "time_list = first_col_dict['Time_(sec)']\n",
    "Plist = first_col_dict['Pressure_(bar)']\n",
    "Tlist = first_col_dict['Temperature_(k)']\n",
    "Vlist = first_col_dict['Volume_(cm3)']\n",
    "\n",
    "print(\"Investigated time point is {0} secs\".format(time_list[timepoint_index]))\n",
    "print(\"At this moment, the reactant's mole fraction is {}.\".format(spc_mf_dict[reactant_name][timepoint_index]))\n",
    "\n",
    "# Create mole_fraction_dict at the investigated time point\n",
    "mole_fraction_dict = {}\n",
    "for i, v in spc_mf_dict.items():\n",
    "    mole_fraction_dict.update({i: v[timepoint_index]})\n",
    "\n",
    "# Create species_identifier_didct for the core species    \n",
    "species_identifier_dict = {}\n",
    "for i, s in enumerate(model.core.species):\n",
    "    species_identifier_dict.update({get_species_identifier(model.core.species[i]): s})\n",
    "\n",
    "# find the total consumption of the reactant in mol (positive value means consumption)\n",
    "reactant_consumption = 0 # \n",
    "for t in range(timepoint_index):\n",
    "    reactantROP0 = spc_total_rop_dict[reactant_name][1][t] # in mole/cm3*s\n",
    "    reactantROP1 = spc_total_rop_dict[reactant_name][1][t+1] # in mole/cm3*s\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "    reactant_consumption += 0.5 * (reactantROP1 * Vlist[t+1] + reactantROP0 * Vlist[t]) * (time_list[t+1] - time_list[t])\n",
    "reactant_consumption = reactant_consumption * -1.0 # in mole\n",
    "\n",
    "if show_integrated_flux:\n",
    "    print(\"At this moment, the reactant's total molar consumption is {0:.3E} mole.\".format(reactant_consumption))\n",
    "\n",
    "# create DiGraph for this moment\n",
    "G = nx.DiGraph()\n",
    "for pair in graph_edges_dict:\n",
    "    node1 = get_species_identifier(pair[0])\n",
    "    node2 = get_species_identifier(pair[1])\n",
    "    e_rawdata = graph_edges_dict[pair]\n",
    "    total_flux = 0\n",
    "    for rxn in e_rawdata:\n",
    "        if show_integrated_flux:\n",
    "            for t in range(timepoint_index):\n",
    "                rxnROP0 = e_rawdata[rxn][t] # in mole/cm3\n",
    "                rxnROP1 = e_rawdata[rxn][t+1] # in mole/cm3\n",
    "                #integrate the flux using a trapezoidal rule\n",
    "                total_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "        else:\n",
    "            total_flux += e_rawdata[rxn][timepoint_index] # in mole/cm3*s\n",
    "    if total_flux >= 0:\n",
    "        G.add_edge(node2, node1, **{\"total_flux\": total_flux}) # in G, positive means production of node1\n",
    "    else:\n",
    "        G.add_edge(node1, node2, **{\"total_flux\": -total_flux}) # in G, negative means consumption of node1   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 1: find most dominant species at the investigated time point\n",
    "\n",
    "Pleae specify\n",
    "\n",
    "- dominant_list: If you want to visualize the first top 10 dominant species, specify the list as [0, 9]. If you want to visualize from the 5th most dominant species to the 10th most dominant species, specify the list as [4, 9].\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "dominant_list = [0, 9] \n",
    "##########################################\n",
    "\n",
    "\n",
    "# sort out the mole_fraction_dict in the descending order\n",
    "v = list(mole_fraction_dict.values())\n",
    "k = list(mole_fraction_dict.keys())\n",
    "index_list = [i[0] for i in sorted(enumerate(v), key=lambda x:x[1], reverse=True)]\n",
    "print('==============================================================')\n",
    "print('Order \\t Species \\t Mole Frac. \\t MW (g/mol)')\n",
    "print('==============================================================')\n",
    "\n",
    "for order in range(dominant_list[0], dominant_list[-1]+1):\n",
    "    index_number = index_list[order]\n",
    "    species_dominant = k[index_number]\n",
    "    molecule_dominant = species_identifier_dict.get(species_dominant)\n",
    "    mole_fraction = v[index_number]\n",
    "    molecular_weight = molecule_dominant.molecule[0].get_molecular_weight()*1000\n",
    "    if len(species_dominant) < 6: # if the species label is too short, put more space to allign it\n",
    "        print('{0} \\t {1} \\t\\t {2} \\t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight))\n",
    "    else:\n",
    "        print('{0} \\t {1} \\t {2} \\t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight))\n",
    "    display(molecule_dominant)\n",
    "    print('--------------------------------------------------------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 2: find paths between two species\n",
    "\n",
    "Pleae specify\n",
    "\n",
    "- source species\n",
    "\n",
    "- target species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "src_species = 'ethane(1)'\n",
    "tgt_species = 'C#C(25)'\n",
    "##########################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## output1: sort out all the paths and most significant one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paths = list(nx.all_simple_paths(G, source=src_species, target=tgt_species, cutoff=5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_fluxes = []\n",
    "if show_integrated_flux:\n",
    "    print('========================================================================================================')\n",
    "    print('Path index \\t Integrated Flux (mole)    Path steps')\n",
    "else:\n",
    "    print('================================================')\n",
    "    print('Path index \\t Flux (mole/cm3/s)    Path steps')\n",
    "print('')\n",
    "\n",
    "for i, path in enumerate(paths):\n",
    "    \n",
    "    path_steps = len(path) - 1\n",
    "    fluxes = [G[path[step]][path[step+1]]['total_flux'] for step in range(path_steps) ]\n",
    "    path_fluxes.append(min(fluxes))\n",
    "    if show_integrated_flux:\n",
    "        print('Path #{0}: \\t {1:.3E} \\t\\t  {2}'.format(i, min(fluxes), path))\n",
    "    else:\n",
    "        print('Path #{0}: \\t {1:.3E} \\t {2}'.format(i, min(fluxes), path))\n",
    "sorted_path_fluxes = sorted(path_fluxes)\n",
    "\n",
    "print('')\n",
    "print('=======================================================')\n",
    "print('\\t Most Significant Path ')\n",
    "print('=======================================================')\n",
    "if show_integrated_flux:\n",
    "    print('Path index \\t Integrated Flux (mole)')\n",
    "    print('Path #{0}: \\t {1:.3E}'.format(path_fluxes.index(sorted_path_fluxes[-1]), sorted_path_fluxes[-1]))\n",
    "else:\n",
    "    print('Path index \\t Flux (mole/cm3/s)')\n",
    "    print('Path #{0}: \\t {1:.3E} '.format(path_fluxes.index(sorted_path_fluxes[-1]), sorted_path_fluxes[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [user input] output2: visualize the path of interest\n",
    "\n",
    "Please specify\n",
    "\n",
    "- the path index you want to visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "path_index_investigate = 0\n",
    "##########################################\n",
    "\n",
    "path = paths[path_index_investigate]\n",
    "path_steps = len(path) - 1\n",
    "print('')\n",
    "print('\\t Pathway Report ')\n",
    "print('======================================')\n",
    "print('The pathway you are intested in has {0} steps.'.format(len(path)))\n",
    "for step in range(path_steps):\n",
    "    step_pair = (path[step], path[step+1])\n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    \n",
    "    print(\"\")\n",
    "    if show_integrated_flux:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1))\n",
    "        print('***********************************************************************')\n",
    "        print(\"{1} --> {2} \\t\\t {3:.3E}\"\n",
    "              .format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux']))\n",
    "        print('***********************************************************************')\n",
    "    else:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1))\n",
    "        print('*****************************************************')\n",
    "        print(\"{1} --> {2} \\t\\t {3:.3E}\"\n",
    "              .format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux']))\n",
    "        print('*****************************************************')\n",
    "    \n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (step_pair[1], step_pair[0])\n",
    "    \n",
    "    print('')\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print(\"Example reaction: rxn#{0}\".format(rxn.index))\n",
    "            print('')\n",
    "            print(str(rxn))\n",
    "            display(rxn)\n",
    "            print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.get_rate_coefficient(T)))\n",
    "            reverse_rate = rxn.generate_reverse_rate_coefficient()\n",
    "            print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "            print('')\n",
    "    if len(h_abs_rxns) > 0: \n",
    "        \n",
    "        print(\"Example H_Abstraction: rxn#{0} (1/{1} H_Abs): \".format(h_abs_rxns[0].index, len(h_abs_rxns)) )\n",
    "        print('')\n",
    "        print(str(h_abs_rxns[0]))\n",
    "        display(h_abs_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = h_abs_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(disp_rxns) > 0: \n",
    "        \n",
    "        print(\"Example Disproportionation: rxn#{0} (1/{1} Disp): \".format(disp_rxns[0].index, len(disp_rxns)) )\n",
    "        print('')\n",
    "        print(str(disp_rxns[0]))\n",
    "        display(disp_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = disp_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        \n",
    "        print(\"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) )\n",
    "        print('')\n",
    "        print(str(Pdep_rxns[0]))\n",
    "        display(Pdep_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].get_rate_coefficient(T, P)))\n",
    "        reverse_rate = Pdep_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverse_rate.get_rate_coefficient(T, P)))\n",
    "        print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 3: find paths from a particular species\n",
    "\n",
    "Please specify\n",
    "\n",
    "- source species\n",
    "\n",
    "- depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps down from source species\n",
    "\n",
    "- path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "source = \"ethane(1)\"\n",
    "depth = 2\n",
    "path_top_list = [0, 0]\n",
    "##########################################\n",
    "\n",
    "\n",
    "current_node = source\n",
    "\n",
    "print('')\n",
    "print('\\t Pathway Report ')\n",
    "print('======================================')\n",
    "print('The pathway you are intested in has {0} steps.'.format(depth))\n",
    "for step in range(depth):\n",
    "    print(\"\\n\")\n",
    "    next_node_flux_list = [(next_node, G[current_node][next_node]['total_flux']) for next_node in G[current_node]]\n",
    "    sorted_next_node_flux_list = sorted(next_node_flux_list, key=lambda tup: -tup[1])\n",
    "    \n",
    "    # choose the top one as next node\n",
    "    tup = sorted_next_node_flux_list[path_top_list[step]]\n",
    "    next_node = tup[0]\n",
    "    step_flux = tup[1]\n",
    "    \n",
    "    print(\"\")\n",
    "    if show_integrated_flux:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1))\n",
    "        print('************************************************************************')\n",
    "        print(\"{0} --> {1} \\t\\t {2:.3E}\"\n",
    "              .format(current_node, next_node, step_flux))\n",
    "        print('************************************************************************')\n",
    "    else:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1))\n",
    "        print('****************************************************')\n",
    "        print(\"{0} --> {1} \\t\\t {2:.3E}\"\n",
    "              .format(current_node, next_node, step_flux))\n",
    "        print('****************************************************')\n",
    "\n",
    "    \n",
    "    step_pair = (current_node, next_node)\n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (next_node, current_node)\n",
    "    \n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print('')\n",
    "            print(\"Example reaction: rxn#{0}\".format(rxn.index))\n",
    "            print('')\n",
    "            if show_integrated_flux:\n",
    "            #integrate the flux using a trapezoidal rule\n",
    "                integrated_flux = 0\n",
    "                for t in range(timepoint_index):\n",
    "                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "                print(\"{0}: Integrated Flux = {1:.3E} mole\".format(str(rxn), integrated_flux))\n",
    "            else:\n",
    "                print(\"{0}: Flux = {1:.3E} mole/cm3/s\".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index]))\n",
    "            display(rxn)\n",
    "            print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.get_rate_coefficient(T)))\n",
    "            reverse_rate = rxn.generate_reverse_rate_coefficient()\n",
    "            print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "            print('')\n",
    "    if len(h_abs_rxns) > 0:\n",
    "        print('')\n",
    "        print(\"Example H_Abstraction: rxn#{0}(1/{1} H_Abs)\".format(h_abs_rxns[0].index, len(h_abs_rxns)))\n",
    "        print('')\n",
    "        print(str(h_abs_rxns[0]))\n",
    "        display(h_abs_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = h_abs_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(disp_rxns) > 0: \n",
    "        print('')\n",
    "        print(\"Example Disproportionation: rxn#{0}(1/{1} Disp)\".format(disp_rxns[0].index, len(disp_rxns)))\n",
    "        print('')\n",
    "        print(str(disp_rxns[0]))\n",
    "        display(disp_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = disp_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        print('')\n",
    "        print(\"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) )\n",
    "        print('')\n",
    "        print(str(Pdep_rxns[0]))\n",
    "        display(Pdep_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].get_rate_coefficient(T, P)))\n",
    "        reverse_rate = Pdep_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverse_rate.get_rate_coefficient(T, P)))\n",
    "        print('')\n",
    "    \n",
    "    current_node = next_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More detail for the H_Abstraction and Disproportionation in last step of the pathway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('')\n",
    "print('The step you see currently is between the two species {0}'.format(step_pair))\n",
    "print('')\n",
    "print('Positive flux in below means positive production of 1st node of the pair.')\n",
    "print('')\n",
    "print('This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns)))\n",
    "\n",
    "flux_threshold = 1e-9\n",
    "\n",
    "print('')\n",
    "print('TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold))\n",
    "\n",
    "total_flux = 0\n",
    "rxn_flux_tups = []\n",
    "for rxn in h_abs_rxns + disp_rxns:\n",
    "    if show_integrated_flux:\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "        flux = 0\n",
    "        for t in range(timepoint_index):\n",
    "            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "    else:\n",
    "        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]\n",
    "    rxn_flux_tups.append((rxn, flux))\n",
    "    \n",
    "\n",
    "rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)\n",
    "for tup in rxn_flux_tups:\n",
    "    rxn = tup[0]\n",
    "    flux = tup[1]\n",
    "    if flux > flux_threshold:\n",
    "        total_flux += flux\n",
    "        print('')\n",
    "        print(\"**********************************************************************************\")\n",
    "        if show_integrated_flux:\n",
    "            print(\"rxn#{0}: {1}: Integrated Flux = {2:.3E} mole \".format(rxn.index, str(rxn), flux) )\n",
    "        else:\n",
    "            print(\"rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s\".format(rxn.index, str(rxn), flux))\n",
    "        display(rxn) \n",
    "print(\"***************************************\")\n",
    "print('')\n",
    "if show_integrated_flux:\n",
    "    print(\"TOTAL integrated flux from h_abs and disp is {0:.3E} mole.\".format(total_flux))\n",
    "else:\n",
    "    print(\"TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.\".format(total_flux))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 4: find paths towards a particular species\n",
    "\n",
    "Similar to functionality 2, please specify\n",
    "\n",
    "Please specify\n",
    "\n",
    "- target species\n",
    "\n",
    "- depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps up beyond target species\n",
    "\n",
    "- path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "target = \"C#C(25)\"\n",
    "depth = 2\n",
    "path_top_list = [0, 0]\n",
    "##########################################\n",
    "\n",
    "\n",
    "current_node = target\n",
    "print('')\n",
    "print('\\t Pathway Report ')\n",
    "print('======================================')\n",
    "print('The pathway you are intested in has {0} steps.'.format(depth))\n",
    "print('')\n",
    "if show_integrated_flux:\n",
    "    product_flux = 0\n",
    "    for t in range(timepoint_index):\n",
    "        productROP0 = spc_total_rop_dict[target][1][t]\n",
    "        productROP1 = spc_total_rop_dict[target][1][t+1]\n",
    "        product_flux += 0.5 * (productROP1 * Vlist[t+1] + productROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "    print(\"total integrated product flux for {0} is {1:.3E} mole.\".format(target, product_flux))\n",
    "else:\n",
    "    print(\"total product flux for {0} is {1:.3E} mole/cm3/s.\".format(target, spc_total_rop_dict[target][1][timepoint_index]))\n",
    "for step in range(depth):\n",
    "    print(\"\\n\")\n",
    "    prev_nodes = []\n",
    "    for node1 in G:\n",
    "        if current_node in G[node1]:\n",
    "            prev_nodes.append(node1)\n",
    "    prevNode_flux_list = [(prev_node, G[prev_node][current_node]['total_flux']) for prev_node in prev_nodes]\n",
    "    sorted_prevNode_flux_list = sorted(prevNode_flux_list, key=lambda tup: -tup[1])\n",
    "    \n",
    "    # choose the top one as next node\n",
    "    tup = sorted_prevNode_flux_list[path_top_list[step]]\n",
    "    prev_node = tup[0]\n",
    "    step_flux = tup[1]\n",
    "    \n",
    "    print(\"\")\n",
    "    if show_integrated_flux:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1))\n",
    "        print('************************************************************************')\n",
    "        print(\"{0} <-- {1} \\t\\t {2:.3E}\"\n",
    "              .format(current_node, next_node, step_flux))\n",
    "        print('************************************************************************')\n",
    "    else:\n",
    "        print(\"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1))\n",
    "        print('****************************************************')\n",
    "        print(\"{0} <-- {1} \\t\\t {2:.3E}\"\n",
    "              .format(current_node, next_node, step_flux))\n",
    "        print('****************************************************')\n",
    "    \n",
    "    step_pair = (prev_node, current_node)\n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (current_node, prev_node)\n",
    "    \n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print('')\n",
    "            print(\"Example reaction: rxn#{0}\".format(rxn.index))\n",
    "            print('')\n",
    "            if show_integrated_flux:\n",
    "            #integrate the flux using a trapezoidal rule\n",
    "                integrated_flux = 0\n",
    "                for t in range(timepoint_index):\n",
    "                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "                print(\"{0}: Integrated Flux = {1:.3E} mole\".format(str(rxn), integrated_flux))\n",
    "            else:\n",
    "                print(\"{0}: Flux = {1:.3E} mole/cm3/s\".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index]))\n",
    "            display(rxn)\n",
    "            print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.get_rate_coefficient(T)))\n",
    "            reverse_rate = rxn.generate_reverse_rate_coefficient()\n",
    "            print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "            \n",
    "    if len(h_abs_rxns) > 0:\n",
    "        print('')\n",
    "        print(\"Example H_Abstraction: rxn#{0}(1/{1} H_Abs)\".format(h_abs_rxns[0].index, len(h_abs_rxns)))\n",
    "        print('')\n",
    "        print(str(h_abs_rxns[0]))\n",
    "        display(h_abs_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = h_abs_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(disp_rxns) > 0: \n",
    "        print('')\n",
    "        print(\"Example Disproportionation: rxn#{0}(1/{1} Disp)\".format(disp_rxns[0].index, len(disp_rxns)))\n",
    "        print('')\n",
    "        print(str(disp_rxns[0]))\n",
    "        display(disp_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].get_rate_coefficient(T)))\n",
    "        reverse_rate = disp_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverse_rate.get_rate_coefficient(T)))\n",
    "        print('')\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        print('')\n",
    "        print(\"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) )\n",
    "        print('')\n",
    "        print(str(Pdep_rxns[0]))\n",
    "        display(Pdep_rxns[0])\n",
    "        print('Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].get_rate_coefficient(T, P)))\n",
    "        reverse_rate = Pdep_rxns[0].generate_reverse_rate_coefficient()\n",
    "        print('Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverse_rate.get_rate_coefficient(T, P)))\n",
    "        print('')\n",
    "    \n",
    "    current_node = prev_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More detail for the H_Abstration and Disproportionation in last step of the pathway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('')\n",
    "print('The step you see currently is between the two species {0}'.format(step_pair))\n",
    "print('')\n",
    "print('Positive flux in below means positive production of 1st node of the pair.')\n",
    "print('')\n",
    "print('This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns)))\n",
    "\n",
    "flux_threshold = 1e-9\n",
    "\n",
    "print('')\n",
    "print('TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold))\n",
    "\n",
    "total_flux = 0\n",
    "rxn_flux_tups = []\n",
    "for rxn in h_abs_rxns + disp_rxns:\n",
    "    if show_integrated_flux:\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "        flux = 0\n",
    "        for t in range(timepoint_index):\n",
    "            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (time_list[t+1] - time_list[t]) # in mole\n",
    "    else:\n",
    "        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]\n",
    "    rxn_flux_tups.append((rxn, flux))\n",
    "\n",
    "rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)\n",
    "for tup in rxn_flux_tups:\n",
    "    rxn = tup[0]\n",
    "    flux = tup[1]\n",
    "    if flux > flux_threshold:\n",
    "        total_flux += flux\n",
    "        print('')\n",
    "        print(\"**********************************************************************************\")\n",
    "        if show_integrated_flux:\n",
    "            print(\"rxn#{0}: {1}: Integrated Flux = {2:.3E} mole \".format(rxn.index, str(rxn), flux) )\n",
    "        else:\n",
    "            print(\"rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s\".format(rxn.index, str(rxn), flux))\n",
    "        display(rxn) \n",
    "print(\"***************************************\")\n",
    "print('')\n",
    "if show_integrated_flux:\n",
    "    print(\"TOTAL integrated flux from h_abs and disp is {0:.3E} mole.\".format(total_flux))\n",
    "else:\n",
    "    print(\"TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.\".format(total_flux))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:rmg_env]",
   "language": "python",
   "name": "conda-env-rmg_env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
